5 de Mayo de 2022
Entender cómo usar los datos geolocalizados para problemas de marketing
En especial
Explicar la demografía de los visitantes de cada centro comercial
Usar los datos de movimientos de la gente para caracterizar la atracción de cada centro comercial
Estudiar problemas de logística/cierre centros comerciales
Recordamos: geomarketing es un tipo de “location intelligence” que, utilizando datos geolocalizados de tiendas, clientes y competencia permite:
El objetivo de este Taller es encontrar las zonas de atracción en los barrios de Madrid para cada uno de los Ikeas (u otros centros comerciales)
Para ello:
Vamos a utilizar datos geolocalizados de Twitter para determinar donde viven los usuarios que van a cada uno de los centros de Ikea de Madrid
Hacemos un modelo (fit) para determinar mediante el modelo de Huff cual es la probabilidad \(p_{ij}\) de que un usuario en un barrio \(i\) visite el centro de Ikea \(j\)
Con ese modelo determinamos las zonas de atracción de cada Ikea, asignando cada barrio \(i\) a un Ikea \(j\) si la probabilidad \(p_{ij}\) es mayor de las \(p_{i*}\)
Vamos a utilizar una base de datos de tweets geolocalizados durante 2 años
load("./data/tweets_madrid.RData")
head(tweets_madrid)
## lat lon user time ## 21 40.49162 -3.689221 1 2012-11-30 19:03:40 ## 22 40.41296 -3.723460 2 2012-11-30 19:05:52 ## 36 40.41916 -3.706140 3 2012-11-30 19:03:07 ## 48 40.39049 -3.725962 4 2012-11-30 19:04:01 ## 54 40.42712 -3.673200 5 2012-11-30 19:07:19 ## 63 40.40596 -3.683690 6 2012-11-30 19:04:22
Cargamos los shapefiles de los barrios de Madrid
map_barrios <- read_sf("./data/Barrios Madrid/200001469.shp")
map_barrios$DESBDT <- iconv(map_barrios$DESBDT,from="Windows-1252", to="UTF-8")
Los transformamos a la proyección WSG84
mapb <- st_transform(map_barrios,4326)
Los visualizamos
leaflet(mapb) %>% addTiles() %>% addPolygons(weight=1)
También necesitaremos la población de los barrios de Madrid
http://www-2.munimadrid.es/TSE6/control/seleccionDatosBarrio
pob_barrios <- read.csv("./data/Madrid-Poblacion_Barrios.csv",
colClasses=c("character","character","integer",
"character","character"))
Añadimos a la tabla del mapa una columna que sea la población (uniendo por el GEOCODIGO)
mapb <- merge(mapb,pob_barrios,by="GEOCODIGO") head(mapb)
## Simple feature collection with 6 features and 7 fields ## Geometry type: POLYGON ## Dimension: XY ## Bounding box: xmin: -3.722965 ymin: 40.40503 xmax: -3.690376 ymax: 40.43061 ## Geodetic CRS: WGS 84 ## GEOCODIGO CODBDT DESBDT Distrito Barrio Población Fecha ## 1 079011 797664 011 Palacio CENTRO PALACIO 22471 10/1/2013 ## 2 079012 797665 012 Embajadores CENTRO EMBAJADORES 46264 10/1/2013 ## 3 079013 797666 013 Cortes CENTRO CORTES 10717 10/1/2013 ## 4 079014 797667 014 Justicia CENTRO JUSTICIA 16847 10/1/2013 ## 5 079015 797668 015 Universidad CENTRO UNIVERSIDAD 31714 10/1/2013 ## 6 079016 797669 016 Sol CENTRO SOL 7577 10/1/2013 ## geometry ## 1 POLYGON ((-3.709384 40.4224... ## 2 POLYGON ((-3.702836 40.4139... ## 3 POLYGON ((-3.696961 40.4189... ## 4 POLYGON ((-3.699416 40.4284... ## 5 POLYGON ((-3.712235 40.4302... ## 6 POLYGON ((-3.70416 40.42019...
En concreto los centros comerciales en la comunidad de Madrid (http://www.madrid.org/nomecalles/)
centros <- read.table("./data/centros_comerciales.csv",header=T,stringsAsFactors = F)
head(centros,2)
## CMUN ETIQUETA MUNICIPIO ## 1 5 Leroy Merlin Alcalá de Henares ## 2 7 Leroy Merlin Alcorcón ## DIRECCIÓN ## 1 Ctra Nacional II, km. 34 - Calle Tamajón s/n (Centro Comercial La Dehesa) ## 2 Avenida de Europa, 26 (Centro Comercial Parque Oeste - Ctra. Nacional V, km. 14) ## lon lat ## 1 -3.327188 40.50292 ## 2 -3.849481 40.34330
Estimamos donde está el domicilio de cada usuario a nivel de barrio: barrio más frecuentado para hacer los tweets.
Para ello añadimos una columna GEOCODIGO en la tabla de tweets de cada barrio:
points_tweets <- st_as_sf(tweets_madrid,coords=c("lon","lat"),crs=4326)
id_polygon <- st_within(points_tweets,mapb)
tweets_madrid$GEOCODIGO <- mapb$GEOCODIGO[as.numeric(id_polygon)]
head(tweets_madrid)
Si va muy lento podemos utilizar también el paquete sp
WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
points_tweets <- cbind(tweets_madrid$lon,tweets_madrid$lat)
pps1 <- SpatialPoints(points_tweets,proj4str=WGS84)
oo1 <- over(pps1,as(as_Spatial(mapb),"SpatialPolygons"))
tweets_madrid$GEOCODIGO <- mapb$GEOCODIGO[oo1]
head(tweets_madrid)
O podemos utilizar los tweets que ya están procesados en
load("./data/tweets_madrid_geo.RData")
Nos quedamos solo con los tweets que suceden en alguno de los barrios
tws <- tweets_madrid[!is.na(tweets_madrid$GEOCODIGO),]
Es posible que haya usuarios que han venido a Madrid de vacaciones y hayan hecho muchos tweets. Nos quedamos solo con los tweets de aquellos usuarios que al menos haya estado en Madrid más de 5 meses. También nos quedamos sólo con aquellos que hayan hecho más de 10 tweets y menos de 10000.
tws$month <- month(tws$time)
users_table <- tws %>%
group_by(user) %>%
summarize(ntot=length(user),nmonths=length(unique(month)))
users_sel <- users_table$user[users_table$nmonths > 5 & users_table$ntot > 10 & users_table$ntot < 10000]
tws <- tws[tws$user %in% users_sel,]
Finalmente calculamos cual ha sido el barrio más frecuentado y cuánto ha sido frecuentado. Para cada usuario obtenemos cuantos tweets tiene en cada barrio
user_barrio_tweets <- tws %>% group_by(user,GEOCODIGO) %>% summarize(nt = length(user))
Ordenamos y calculamos el barrio más visitado
user_barrio_tweets <- user_barrio_tweets[order(user_barrio_tweets$user,user_barrio_tweets$nt,decreasing=T),]
user_barrio_most <- user_barrio_tweets %>%
group_by(user) %>%
summarize(ntot = sum(nt),
nmax=head(nt,1),
bmax=head(GEOCODIGO,1),
nbarrios=length(unique(GEOCODIGO)))
Solo consideramos aquellos en los que más del 30% de los tweets suceden en el barrio más visitado
user_barrio_most <- filter(user_barrio_most,ntot < 10000 & ntot > 10 & nmax/ntot > 0.3 & nbarrios > 3)
Vamos a ver cómo se comparan el número de homes con la población de cada barrio
aa <- user_barrio_most %>% group_by(bmax) %>% summarise(nusers=length(user)) %>%
ungroup()
merge(aa,pob_barrios,by.x="bmax",by.y="GEOCODIGO") %>%
ggplot(aes(x=nusers,y=Población)) + geom_point() +
xlab("Twitter users") + ylab("Población")
Estimamos cuantos usuarios han ido a un centro comercial, desde que barrio y a qué distancia.
Vamos a probar con los de Ikea. Con el de Vallecas
centros_sel <- centros[centros$ETIQUETA=="Ikea",]
points_IKEA <- st_as_sf(centros_sel,coords=c("lon","lat"),crs=4326)
Vemos qué usuarios han geolocalizado un tweet a menos de 500 metros del centro de IKEA
points_tweets <- st_as_sf(tws,coords=c("lon","lat"),crs=4326)
distancia_a_Ikea <- st_distance(points_tweets,points_IKEA)
users_Ikea <- unique(tws$user[as.numeric(distancia_a_Ikea[,3]) < 500])
Hacemos una tabla de los barrios de donde vienen esos usuarios y la distancia (al centro del barrio)
users_barrio_Ikea <- user_barrio_most[user_barrio_most$user %in% users_Ikea,] agg_barrio_Ikea <- users_barrio_Ikea %>% group_by(bmax) %>% summarise(nn=length(user))
Juntamos estos datos con el mapa
mapb_Ikea <- mapb mapb_Ikea <- merge(mapb_Ikea,agg_barrio_Ikea,by.x="GEOCODIGO",by.y="bmax")
Pintamos un choropleth con estos datos
mapp <- mapb_Ikea[!is.na(mapb_Ikea$nn),]
pal <- colorNumeric("magma",domain=mapp$nn)
m <- leaflet(mapp) %>% addTiles() %>%
addPolygons(fillColor= ~ pal(nn),
weight=1,fillOpacity = 0.8,
popup=~paste(DESBDT,"<br> nusers = ",nn,sep=""))
Veamos como es este choropleth
m
Ejercicio 1:
¿Desde qué barrio van más personas al centro de Ikea del norte?
¿Cuál es la distancia media que recorren los clientes para ir a cada uno de los Ikeas?
Finalmente, vamos a modelizar la atracción a cada uno de los centros comerciales (Ikea). Para ello vamos a utilizar la ley de Huff que dice que la probabilidad de que una persona que vive en \(i\) visite un centro comercial en \(j\) viene dada por
\[ p_{ij} \sim \frac{A_j}{d_{ij}^\beta} \]
donde \(A_j\) es lo atractivo que es el centro comercial \(j\) (normalmente se toma el tamaño), mientras que \(d_{ij}\) es la distancia del cliente \(i\) al centro \(j\).
En nuestro caso vamos a empezar por ver cómo la distancia determina el número de usuarios que va a un centro comercial, mediante la ley
\[N_{ij} \sim \frac{1}{d_{ij}^\beta}\]
Para ello calculamos primero la distancia del Ikea al centro de cada barrio
cc <- st_centroid(mapp) dd <- as.numeric(st_distance(cc,points_IKEA)[,3])
Hacemos el fit (en logaritmos ya que es una función no lineal)
fit <- lm(log(nn) ~ log(dd), data=mapp) summary(fit)
## ## Call: ## lm(formula = log(nn) ~ log(dd), data = mapp) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1.54144 -0.53346 -0.01751 0.50824 1.61527 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 15.4195 1.5635 9.862 < 2e-16 *** ## log(dd) -1.5286 0.1705 -8.963 1.21e-14 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.6921 on 106 degrees of freedom ## Multiple R-squared: 0.4311, Adjusted R-squared: 0.4258 ## F-statistic: 80.33 on 1 and 106 DF, p-value: 1.207e-14
donde vemos que \(\beta \approx\) 1.53
Pero el numero de personas también depende de la población del barrio, asi que vamos a suponer una ley de la forma
\[N_{ij} \sim \frac{P_i^\alpha}{d_{ij}^\beta}\]
con lo que el modelo sería
fit1 <- lm(log(nn) ~ log(Población) + log(dd), data=mapp) summary(fit1) ## ## Call: ## lm(formula = log(nn) ~ log(Población) + log(dd), data = mapp) ## ## Residuals: ## Min 1Q Median 3Q Max ## -1.48464 -0.48147 0.01948 0.45129 1.65638 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 11.35004 1.77123 6.408 4.27e-09 *** ## log(Población) 0.39113 0.09633 4.060 9.46e-05 *** ## log(dd) -1.51331 0.15935 -9.497 8.26e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 0.6465 on 105 degrees of freedom ## Multiple R-squared: 0.5083, Adjusted R-squared: 0.499 ## F-statistic: 54.28 on 2 and 105 DF, p-value: < 2.2e-16
Dado que el fit es a un número de usuarios (variable discreta), probamos un GLM poissoniano con link logaritmico
\[\log (\lambda_{ij}) \sim \alpha log(P_i) - \beta log(d_{ij})\]
fit2 <- glm(nn ~ log(Población) + log(dd),data=mapp,family=poisson(link="log")) summary(fit2) ## ## Call: ## glm(formula = nn ~ log(Población) + log(dd), family = poisson(link = "log"), ## data = mapp) ## ## Deviance Residuals: ## Min 1Q Median 3Q Max ## -5.1388 -1.1508 -0.4280 0.5194 4.7271 ## ## Coefficients: ## Estimate Std. Error z value Pr(>|z|) ## (Intercept) 8.85963 1.01115 8.762 <2e-16 *** ## log(Población) 0.76976 0.07101 10.841 <2e-16 *** ## log(dd) -1.64339 0.06243 -26.324 <2e-16 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## (Dispersion parameter for poisson family taken to be 1) ## ## Null deviance: 1065.8 on 107 degrees of freedom ## Residual deviance: 246.5 on 105 degrees of freedom ## AIC: 610.31 ## ## Number of Fisher Scoring iterations: 5
Los resultados anteriores nos dan una forma funcional para el número de usuarios del barrio \(i\) que visitan el centro comercial \(j\):
\[N_{ij} \sim \frac{P_i^\alpha}{d_{ij}^\beta}\].
Vamos a utilizarla para determinar las zonas de atracción suponiendo que esa ley vale para todos los centros comerciales
nij_Huff_model <- function(pts1,pts2,pob,alpha,beta){
dd <- st_distance(pts1,pts2)
prob_Huff_model <- pob^alpha/as.numeric(dd)^beta
prob_Huff_model
}
coord_barrios <- st_centroid(mapb)
nij <- data.frame()
for(i in 1:nrow(points_IKEA)){
nij <- rbind(nij,nij_Huff_model(coord_barrios,points_IKEA[i,],
mapb$Población,0.31,1.44))
}
colnames(nij) <- mapb$GEOCODIGO
head(nij[,1:8])
## 079011 079012 079013 079014 079015 079016 ## 1 2.379124e-05 2.838968e-05 1.680126e-05 1.840339e-05 2.366358e-05 1.570227e-05 ## 2 1.805954e-05 2.234352e-05 1.518181e-05 1.890398e-05 2.249597e-05 1.355679e-05 ## 3 3.284778e-05 4.815319e-05 3.127393e-05 3.351510e-05 3.636477e-05 2.544260e-05 ## 079021 079022 ## 1 2.622335e-05 2.870612e-05 ## 2 1.666106e-05 1.924022e-05 ## 3 3.336490e-05 4.480561e-05
Para cada barrio, nos quedamos con el centro que da más número de visitantes. Eso define la zona de atracción. Y la visualizamos
centro_most <- apply(nij,2,which.max) mapb$centro_most <- centro_most factpal <- colorFactor(topo.colors(3), as.factor(mapb$centro_most)) leaflet(mapb) %>% addTiles() %>% addPolygons(weight=2,color = ~factpal(as.factor(centro_most)))
Repetir el análisis de las zonas de atracción suponiendo que el IKEA de Alcorcón cierra.
Repetir el análisis para los centros de Leroy Merlin